
###########################################################################################\
# Necessary Packages
###########################################################################################\

library(ggplot2); library(data.table); library(dplyr); library(dtplyr); library(zoo); library(stringr); library(tidyr); library(forcats)
is.integer64 = function(x) {class(x) == "integer64"}

###########################################################################################\
# Load in Data 
###########################################################################################\

# The provided Y-14 file is an empty placeholder that replicates the column structure of the actual FR Y-14 data.
# It allows the following code to run without access to the underlying confidential data. 
y14 <- fread("y14_data.csv")

# The provided Y-9C file aggregates term loans and unused commitments across banks in the raw FR Y-9C data. 
y9c_data <- fread("y9c_data.csv")

# This REITs List file replicates the column structure of the file used to identify REITs in the data
# It has no observations but allows the following code to run. 
reits_list = fread("REITS_list.csv") %>%
  mutate(cusip = substr(cusip, 1, 6))

# This SRISK data file replicates the column structure of the SRISK data used in the analysis.
# It has no observations but allows the following code to run. 
srisk = fread("srisk_data.csv") 

# SP 500 Returns for SRISK Analysis
sp500 = fread("sp500_returns.csv") 

# Mapping from NAICS codes to titles
naics_titles <- fread("naics_titles_mapping.csv")

###########################################################################################\
# Important Functions 
###########################################################################################\

# Functions to Clean NAICS Codes/Titles

merge_naics_titles = function(input) {
  input %>%
    left_join(naics_titles,
              by = c("four_digit_naics" = "2012 NAICS US   Code")) %>%
    rename(naics_title = `2012 NAICS US Title`) %>%
    mutate(naics = paste0(naics_title, " (", four_digit_naics, ") ")) %>%
    return()
  
}

clean_industry_codes = function(input) {
  
  input = input %>% lazy_dt() %>% mutate(naics_2012_equivalent = as.character(naics_2012_equivalent)) %>%
    mutate(industry_code_type_cat = as.character(industry_code_type_cat))
  
  input = input %>% mutate(naics_2012_equivalent = ifelse(industry_code_id %in% c(-99, 0, 1, 9999, 99999), "0", naics_2012_equivalent))
  input = input %>% mutate(naics_2012_equivalent = ifelse(is.na(industry_code_id), "0", naics_2012_equivalent))
  
  input = input %>% mutate(industry_code_type_cat = ifelse(nchar(industry_code_id) == 8, "3", industry_code_type_cat))
  
  input = input %>% mutate(industry_code_type_cat = ifelse(nchar(industry_code_id) == 6 & is.na(industry_code_type_cat), "1", industry_code_type_cat))
  
  input = input %>% mutate(industry_code_type_cat = ifelse(nchar(industry_code_id) == 6 & industry_code_type_cat == 2, "1", industry_code_type_cat))
  
  input = input %>% mutate(naics_2012_equivalent = ifelse(! is.na(industry_code_type_cat) & 
                                                            ! is.na(industry_code_id) & 
                                                            industry_code_type_cat == 1 & 
                                                            (! substr(industry_code_id, 1, 2) %in% c("52", "53")), "0", naics_2012_equivalent))
  
  input = input %>% mutate(naics_2012_equivalent = ifelse(is.na(naics_2012_equivalent) & industry_code_type_cat == 2 &
                                                            (! substr(industry_code_id, 1, 2) %in% c("64", "66", "65", "68", "77", "60", "61", "62", "63")), "0", naics_2012_equivalent))
  
  input = input %>% mutate(naics_2012_equivalent = ifelse(is.na(naics_2012_equivalent) & industry_code_type_cat == 3 &
                                                            (! substr(industry_code_id, 1, 2) %in% c("40")), "0", naics_2012_equivalent))
  
  input = input %>% mutate(naics_2012_equivalent = ifelse(is.na(naics_2012_equivalent) & 
                                                            ! is.na(industry_code_id) &
                                                            industry_code_type_cat == 1, substr(industry_code_id, 1, 4), naics_2012_equivalent))
  
  input = input %>% mutate(naics_2012_equivalent = ifelse(industry_code_type_cat == 1 & 
                                                            substr(naics_2012_equivalent, 1, 4) %in% c(5252, 5213), "0", naics_2012_equivalent))
  
  input = input %>% mutate(two_digit_naics = substr(naics_2012_equivalent, 1, 2),
                           four_digit_naics = substr(naics_2012_equivalent, 1, 4)) %>%
    filter(! four_digit_naics %in% c("5221", "5211"))
  
  return(input %>% as.data.frame())
  
}

# Functions for SRISK Granger Casuality Analysis

extract_alpha = function(input) {
  
  bank = input %>% filter(str_detect(series, "_bank_")) %>% mutate(abnormal_return = NA)
  nonbank = input %>% filter(str_detect(series, "_nonbank_")) %>% mutate(abnormal_return = NA)
  
  for (i in 95:NROW(bank)) {
    subset = bank %>% slice((i-94):(i-4))
    mod = lm(index ~ sprtrn, data = subset)
    beta = mod$coefficients["sprtrn"] %>% as.numeric()
    bank$abnormal_return[i] = bank$index[i] - beta*bank$sprtrn[i]
    print(i)
  }
  for (i in 95:NROW(nonbank)) {
    subset = nonbank %>% slice((i-94):(i-4))
    mod = lm(index ~ sprtrn, data = subset)
    beta = mod$coefficients["sprtrn"] %>% as.numeric()
    nonbank$abnormal_return[i] = nonbank$index[i] - beta*nonbank$sprtrn[i]
    print(i)
  }
  
  output = rbind(bank, nonbank)
  return(output)
}



date_to_period_mapping = function(input) {
  input %>%
    mutate(Period = ifelse(date <= "2007-07-31", "Pre-GFC", NA),
           Period = ifelse(date >= "2007-08-01" & date <= "2009-08-31", "GFC", Period),
           Period = ifelse(date >= "2009-09-01" & date <= "2016-06-30", "Post-GFC", Period),
           Period = ifelse(date >= "2014-12-01" & date <= "2016-06-30", "Oil Shock", Period),
           Period = ifelse(date >= "2016-07-01" & date <= "2019-12-31", "Hike + QT", Period),
           Period = ifelse(date >= "2020-01-01" & date <= "2021-10-29", "Pandemic", Period),
           Period = ifelse(date >= "2021-11-01" & date <= "2022-12-30", "Post-Pandemic", Period),
           Period = ifelse(date >= "2023-01-01", "SVB Stress", Period)) %>% 
    mutate(period_number = ifelse(date <= "2007-07-31", 2, NA),
           period_number = ifelse(date >= "2007-08-01" & date <= "2009-08-31", 3, period_number),
           period_number = ifelse(date >= "2009-09-01" & date <= "2016-06-30", 4, period_number),
           period_number = ifelse(date >= "2014-12-01" & date <= "2016-06-30", 5, period_number),
           period_number = ifelse(date >= "2016-07-01" & date <= "2019-12-31", 6, period_number),
           period_number = ifelse(date >= "2020-01-01" & date <= "2021-10-29", 7, period_number),
           period_number = ifelse(date >= "2021-11-01" & date <= "2022-12-30", 8, period_number),
           period_number = ifelse(date >= "2023-01-01", 9, period_number)) %>%
    mutate(Period = factor(Period, levels = c("Pre-GFC", "GFC", "Post-GFC", "Oil Shock", "Hike + QT", "Pandemic", "Post-Pandemic", "SVB Stress"))) %>%
    return()
  
}


granger_function = function(x, y, max_lags, optimal_lags) {
  return(lmtest::grangertest(y ~ x, order = select_lags(x, y, max_lags, optimal_lags))$`Pr(>F)`[2] %>% round(3))
}


select_lags = function(predictor, outcome, max_lags, optimal_lags) {
  
  if (optimal_lags) {
    y = outcome
    x = predictor
    
    out = c()
    for (k in 1:max_lags) {
      
      formula = "y ~ "
      for (i in 1:k) {
        formula = paste0(formula, "lag(y, ", i, ") + ") %>% paste0("lag(x, ", i, ") + ")
      }
      formula = substr(formula, 1, nchar(formula)-2)
      mod = lm(as.formula(formula)) %>% broom::glance()
      
      out = c(out, mod$AIC)
    }
    return(which.min(out))
  }
  else {return(1)}
}


create_correlations = function(input, bank_index, nonbank_index, lag_param, rollcorr_param) {
  
  bank = input %>% filter(series == bank_index) %>% select(date, index) %>% rename(bank = index)
  nonbank = input %>% filter(series == nonbank_index) %>% select(date, index) %>% rename(nonbank = index)
  
  df = inner_join(bank, nonbank, by = "date")
  df = lag_function(df, lag_param, "bank")
  df = df %>% read.zoo() %>% rollapply(width = rollcorr_param, function(x) cor(x[,1], x[,2], use = 'complete.obs'), by.column = FALSE) %>% as.data.frame()
  df$date = rownames(df)
  rownames(df) = NULL
  names(df) = c("cor", "date")
  df$bank = bank_index
  df$nonbank = nonbank_index
  df$Window = paste0(rollcorr_param, " Days")
  
  return(df)
}


lag_function = function(input, param, colname) {
  if (param >= 0) {
    input[[colname]] = lag(input[[colname]], param)
    return(input)
  }
  else {
    input[[colname]] = lead(input[[colname]], param*-1)
    return(input)
  }
}


###########################################################################################\
# Data Cleaning
###########################################################################################\

# Generate a list of banks present for the entire sample (44 quarters as of 2023Q4)
banks_for_y14_balanced_panel = y14 %>% 
  distinct(id_rssd, quarter) %>%
  group_by(id_rssd) %>% mutate(obs = n()) %>% ungroup() %>%
  filter(obs == 44) %>%
  pull(id_rssd)%>% unique()

## Create the underlying data for the Y-14 term loan and credit line graphs
y14_graphing_data = y14 %>% 
  clean_industry_codes() %>%
  filter(two_digit_naics %in% c("52", "53")) %>%
  mutate(four_digit_naics = ifelse(substr(four_digit_naics, 1, 2) == "53", "53", four_digit_naics)) %>%
  mutate(facility_type_cat = as.numeric(facility_type_cat)) %>%
  mutate(facility_type = ifelse(facility_type_cat %in% c(1:6), "Credit Line", ifelse(facility_type_cat %in% c(7:12), "Term Loan", NA))) %>%
  filter(! is.na(facility_type)) %>%
  group_by(quarter, four_digit_naics, facility_type) %>% 
  summarize(committed_exposure_amt = sum(committed_exposure_amt, na.rm = T)) %>% 
  ungroup() %>%
  merge_naics_titles() %>%
  mutate(committed_exposure_amt = committed_exposure_amt/1e9) %>%
  mutate(naics = trimws(naics)) %>%
  filter(naics != "Insurance and Employee Benefit Funds (5251)") %>% filter(naics != "Securities and Commodity Exchanges (5232)")


## Create the underlying data for the graphs related to REITS 
reits_graphing_data = y14 %>%
  lazy_dt() %>%
  filter(id_rssd %in% banks_for_y14_balanced_panel) %>%
  mutate(dt_origination = as.Date(dt_origination, "%d%b%Y")) %>% 
  mutate(origination_quarter = as.yearqtr(dt_origination)) %>%
  mutate(quarter = as.yearqtr(quarter), origination_quarter = as.yearqtr(origination_quarter)) %>%
  filter(committed_exposure_amt > 0) %>%
  filter(facility_type_cat %in% c(1:6)) %>%
  arrange(id_rssd, internal_credit_facility_id, quarter) %>%
  group_by(id_rssd, internal_credit_facility_id) %>%
  mutate(lagged_utilization = lag(utilized_exposure_amt),
         lagged_committed = lag(committed_exposure_amt)) %>%
  ungroup() %>%
  mutate(drawdown = ifelse(origination_quarter < quarter & utilized_exposure_amt >= lagged_utilization, utilized_exposure_amt - lagged_utilization, 0),
         repayment = ifelse(origination_quarter < quarter & utilized_exposure_amt < lagged_utilization, lagged_utilization - utilized_exposure_amt, 0)) %>%
  mutate(counterparty = case_when((tkr %in% reits_list$ticker | issuer_cusip_id %in% reits_list$cusip) ~ "REITs",
                                  (! (tkr %in% reits_list$ticker | issuer_cusip_id %in% reits_list$cusip)) & (substr(naics_2012_equivalent, 1, 2) %in% c("52", "53")) & tkr == "" ~ "Financial Ex-REIT: Private",
                                  (! (tkr %in% reits_list$ticker | issuer_cusip_id %in% reits_list$cusip)) & (! substr(naics_2012_equivalent, 1, 2) %in% c("52", "53")) & tkr == "" ~ "Non-Financial: Private",
                                  (! (tkr %in% reits_list$ticker | issuer_cusip_id %in% reits_list$cusip)) & (substr(naics_2012_equivalent, 1, 2) %in% c("52", "53")) & tkr != "" ~ "Financial Ex-REIT: Public",
                                  (! (tkr %in% reits_list$ticker | issuer_cusip_id %in% reits_list$cusip)) & (! substr(naics_2012_equivalent, 1, 2) %in% c("52", "53")) & tkr != "" ~ "Non-Financial: Public")) %>%
  mutate(`Drawdown Rate` = ifelse(committed_exposure_amt > lagged_committed, drawdown/committed_exposure_amt, drawdown/lagged_committed)) %>%
  mutate(`Net Drawdown Rate` = ifelse(committed_exposure_amt > lagged_committed, (drawdown-repayment)/committed_exposure_amt, (drawdown-repayment)/lagged_committed)) %>%
  mutate(`Utilization Rate` = utilized_exposure_amt/committed_exposure_amt) %>%
  filter(`Utilization Rate` < 1) %>%
  as.data.frame()

# Clean the SRISK data for analysis
srisk_analysis_data = srisk %>% select(date, aggregate_srisk_truncated_bank, aggregate_srisk_truncated_nonbank, 
                                     equity_returns_bank_cap_weighted, equity_returns_nonbank_cap_weighted, 
                                     equity_returns_bank_equal_weighted, equity_returns_nonbank_equal_weighted) %>% 
  pivot_longer(aggregate_srisk_truncated_bank:equity_returns_nonbank_equal_weighted, names_to = "series", values_to = "index") %>% 
  arrange(series, date) %>%
  mutate(log_index = log(index)) %>%
  group_by(series) %>%
  mutate(log_change = log_index - lag(log_index)) %>%
  ungroup()%>%
  mutate(index = ifelse(series %in% c("aggregate_srisk_truncated_bank", "aggregate_srisk_truncated_nonbank"), log_change, index)) %>%
  select(date, series, index)

###########################################################################################\
# Code to Produce Figures in the Paper
###########################################################################################\

############################### FIGURE 2  ###############################

# FIGURE 2A
fig_2a <- y14_graphing_data %>%
  filter(facility_type == "Term Loan") %>%
  mutate(quarter = as.yearqtr(quarter)) %>%
  filter(quarter <= as.yearqtr("2023 Q4")) %>%
  ggplot(aes(x = quarter, y = committed_exposure_amt, fill = naics)) +
  geom_area() + theme_classic() + scale_fill_brewer(palette = "Pastel2") + 
  ggtitle("Term Loan Committed Exposure by NAICS") +
  annotate(geom = "text", x = as.yearqtr("2023 Q3"), y = 60, label = "38%") +
  annotate(geom = "text", x = as.yearqtr("2023 Q3"), y = 190, label = "40%") +
  annotate(geom = "text", x = as.yearqtr("2023 Q3"), y = 280, label = "11%") +
  annotate(geom = "text", x = as.yearqtr("2013 Q2"), y = 25, label = "41%") +
  annotate(geom = "text", x = as.yearqtr("2013 Q2"), y = 68, label = "19%") +
  annotate(geom = "text", x = as.yearqtr("2013 Q2"), y = 92, label = "18%") +
  scale_y_continuous(breaks = c(0, 50, 100, 150, 200, 250, 300)) + xlab("") +
  ylab("Billions of Dollars") + theme(legend.position = "bottom", legend.title = element_blank()) +
  scale_x_yearqtr(format = "%Y Q%q") + guides(fill = guide_legend(nrow = 4, byrow = T))
fig_2a

# FIGURE 2B
fig_2b <- y14_graphing_data %>%
  filter(facility_type == "Credit Line") %>%
  mutate(quarter = as.yearqtr(quarter)) %>%
  filter(quarter <= as.yearqtr("2023 Q4")) %>%
  ggplot(aes(x = quarter, y = committed_exposure_amt, fill = naics)) +
  geom_area() + theme_classic() + scale_fill_brewer(palette = "Pastel2") + 
  ggtitle("Credit Line Committed Exposure by NAICS") +
  annotate(geom = "text", x = as.yearqtr("2023 Q3"), y = 175, label = "20%") +
  annotate(geom = "text", x = as.yearqtr("2023 Q3"), y = 700, label = "43%") +
  annotate(geom = "text", x = as.yearqtr("2023 Q3"), y = 1150, label = "15%")  +
  annotate(geom = "text", x = as.yearqtr("2023 Q3"), y = 1340, label = "11%") +
  annotate(geom = "text", x = as.yearqtr("2013 Q2"), y = 80, label = "24%") +
  annotate(geom = "text", x = as.yearqtr("2013 Q2"), y = 200, label = "28%") +
  annotate(geom = "text", x = as.yearqtr("2013 Q2"), y = 310, label = "13%") +
  annotate(geom = "text", x = as.yearqtr("2013 Q2"), y = 400, label = "19%") +
  scale_y_continuous(breaks = c(0, 500, 1000, 1500)) + xlab("") +
  ylab("Billions of Dollars") + theme(legend.position = "bottom", legend.title = element_blank()) +
  scale_x_yearqtr(format = "%Y Q%q") + guides(fill = guide_legend(nrow = 4, byrow = T))
fig_2b

# FIGURE 2C
fig_2c <-  ggplot(y9c_data, aes(x = quarter, y = share)) +
  geom_line(size = 1.5, color = "steelblue") + 
  facet_wrap(vars(group)) + 
  scale_x_yearqtr(labels = function(x) zoo::format.yearqtr(x, "%Y Q%q")) + 
  xlab("") + ylab("") +
  theme_bw() +
  theme(legend.position = "bottom", legend.title = element_blank()) + 
  scale_y_continuous(labels = scales::percent)
fig_2c

############################### FIGURE 6  ###############################

fig_6 <- reits_graphing_data %>%
  mutate(period = ifelse(quarter >= as.yearqtr("2018 Q1") & quarter <= as.yearqtr("2019 Q4"), "2018:Q1-2019:Q4", ifelse(quarter == as.yearqtr("2020 Q1"), "2020:Q1", NA))) %>%
  mutate(quintile = case_when(pd_dec < 0.0015 ~ "First Quintile: PD < 0.15%",
                              pd_dec >= 0.0015 & pd_dec < 0.0033 ~ "Second Quintile: 0.15% <= PD < 0.33%",
                              pd_dec >= 0.0033 & pd_dec < 0.0067 ~ "Third Quintile: 0.33% <= PD < 0.67%",
                              pd_dec >= 0.0067 & pd_dec < 0.0167 ~ "Fourth Quintile: 0.67% <= PD < 1.67%",
                              pd_dec >= 0.0167 ~ "Fifth Quintile: 1.67% <= PD")) %>%
  mutate(quintile = factor(quintile, levels = c("First Quintile: PD < 0.15%",
                                                "Second Quintile: 0.15% <= PD < 0.33%",
                                                "Third Quintile: 0.33% <= PD < 0.67%",
                                                "Fourth Quintile: 0.67% <= PD < 1.67%",
                                                "Fifth Quintile: 1.67% <= PD"))) %>%
  filter(! is.na(quintile), ! is.na(counterparty), ! is.na(period)) %>%
  group_by(quintile, period, counterparty) %>%
  summarize(`Weighted Average Utilization Rate` = weighted.mean(`Utilization Rate`, w = committed_exposure_amt, na.rm = T) %>% round(3),
            `Weighted Average Drawdown Rate` = weighted.mean(`Drawdown Rate`, w = committed_exposure_amt, na.rm = T) %>% round(3),
            `Weighted Average Net Drawdown Rate` = weighted.mean(`Net Drawdown Rate`, w = committed_exposure_amt, na.rm = T) %>% round(3)) %>%
  ungroup() %>%
  pivot_longer(`Weighted Average Utilization Rate`:`Weighted Average Net Drawdown Rate`, names_to = "measure", values_to = "value") %>%
  pivot_wider(names_from = "counterparty", values_from = "value") %>%
  select(quintile, period, measure, `Financial Ex-REIT: Public`, `Non-Financial: Public`, `REITs`) %>%
  pivot_wider(names_from = "period", values_from = c("Financial Ex-REIT: Public", "Non-Financial: Public", "REITs")) %>%
  arrange(quintile) %>%
  gt() %>%
  tab_spanner(label = "Non-Financial", columns = c("Non-Financial: Public_2018:Q1-2019:Q4", "Non-Financial: Public_2020:Q1")) %>%
  tab_spanner(label = "REITs", columns = c("REITs_2018:Q1-2019:Q4", "REITs_2020:Q1")) %>%
  tab_spanner(label = "Financial Ex-REIT", columns = c("Financial Ex-REIT: Public_2018:Q1-2019:Q4", "Financial Ex-REIT: Public_2020:Q1")) %>%
  cols_label(`Non-Financial: Public_2018:Q1-2019:Q4` = "2018:Q1-2019:Q4",
             `Financial Ex-REIT: Public_2018:Q1-2019:Q4` = "2018:Q1-2019:Q4",
             `REITs_2018:Q1-2019:Q4` = "2018:Q1-2019:Q4",
             `Non-Financial: Public_2020:Q1` = "2020:Q1",
             `Financial Ex-REIT: Public_2020:Q1` = "2020:Q1",
             `REITs_2020:Q1` = "2020:Q1",
             quintile = "Risk Category") %>%
  tab_row_group(label = "Weighted Average Drawdown Rate", rows = measure == "Weighted Average Drawdown Rate") %>%
  tab_row_group(label = "Weighted Average Net Drawdown Rate", rows = measure == "Weighted Average Net Drawdown Rate") %>%
  tab_row_group(label = "Weighted Average Utilization Rate", rows = measure == "Weighted Average Utilization Rate") %>%
  cols_hide(c("measure")) %>%
  row_group_order(c("Weighted Average Utilization Rate", "Weighted Average Drawdown Rate", "Weighted Average Net Drawdown Rate")) %>%
  cols_align(c("center")) %>%
  tab_footnote(footnote = "Notes: Utilization rate defined at the loan level as current period utilized amount divided by current period committed amount. 
  Drawdown rate defined at the loan level as change in utilization divided by previous period committed amount, if the change in utilization is greater than zero. 
  Net drawdown rate defined at the loan level as change in utilization divided by previous period committed amount. Weighted averages are calculated using
  current period committed amount as weights. Default probability quintiles are based on distribution as of 2020:Q1. Underlying data is based on a 
  balanced panel of 23 banks. As of 2020:Q1, underlying data has 1100 unique credit facilities to REITs, 15141 unique credit facilities to public non-financials, 
  66230 unique credit facilities to private non-financials, 2500 unique credit facilities to public non-REIT financials, and 13255 unique credit facilities to private non-REIT financials.")
fig_6

############################### FIGURE 7  ###############################

# FIGURE 7A and 7B (Top Panel)
fig_7a_7b <- reits_graphing_data %>%
  lazy_dt() %>%
  filter(counterparty == "REITs", quarter != as.yearqtr("2013 Q1")) %>%
  group_by(quarter) %>%
  summarize(`Outstanding Commitments` = sum(committed_exposure_amt, na.rm = T),
            `Drawdowns` = sum(drawdown, na.rm = T)) %>%
  ungroup() %>%
  mutate(`Outstanding Commitments` = `Outstanding Commitments`/1e9,
         `Drawdowns` = `Drawdowns`/1e9) %>%
  pivot_longer(`Outstanding Commitments`:`Drawdowns`, names_to = "series", values_to = "values") %>%
  mutate(qtr = as.character(quarter)) %>%
  as.data.frame() %T>%
  mutate(series = factor(series, levels = c("Outstanding Commitments", "Drawdowns"))) %>%
  ggplot(aes(x = quarter, y = values)) +
  geom_line(size = 1.5, color = "steelblue3") + theme_bw() +
  facet_wrap(vars(series), scales = "free_y") +
  xlab("") + ylab("Billions of USD") + scale_x_yearqtr(labels = function(x) zoo::format.yearqtr(x, "%Y Q%q"))
fig_7a_7b

# FIGURE 7C  (Bottom Panel)
fig_7c <-reits_graphing_data %>%
  lazy_dt() %>%
  filter(counterparty %in% c("REITs", "Financial Ex-REIT: Private", "Financial Ex-REIT: Public"), 
         quarter != as.yearqtr("2013 Q1")) %>%
  group_by(counterparty, quarter) %>%
  summarize(committed = sum(committed_exposure_amt, na.rm = T),
            utilized = sum(utilized_exposure_amt, na.rm = T),
            drawdown = sum(drawdown, na.rm = T)) %>%
  ungroup() %>%
  mutate(utilization_rate = utilized/committed) %>%
  select(counterparty, quarter, utilization_rate, drawdown) %>%
  pivot_wider(names_from = "counterparty", values_from = c("utilization_rate", "drawdown")) %>%
  mutate(reit_drawdown_share = drawdown_REITs/(drawdown_REITs + `drawdown_Financial Ex-REIT: Private` + `drawdown_Financial Ex-REIT: Public`)) %>%
  select(quarter, reit_drawdown_share, utilization_rate_REITs) %>%
  as.data.frame() %>%
  setNames(c("quarter", "Drawdown Share of REITs", "REIT Commitment Utilization Rate")) %>%
  pivot_longer(`Drawdown Share of REITs`:`REIT Commitment Utilization Rate`, names_to = "group", values_to = "values") %>%
  mutate(qtr = as.character(quarter)) %>%
  mutate(group = factor(group, levels = c("REIT Commitment Utilization Rate", "Drawdown Share of REITs"))) %>%
  filter(group == "Drawdown Share of REITs") %>%
  ggplot(aes(x = quarter, y = values, values = group)) +
  geom_line(size = 1.5, color = "steelblue3") +
  theme_bw() + xlab("") + ylab("") +
  scale_x_yearqtr(labels = function(x) zoo::format.yearqtr(x, "%Y Q%q")) +
  facet_wrap(vars(group), scales = "free_y") +
  scale_y_continuous(labels = scales::percent)
fig_7c


############################### Figure 9  ###############################

# FIGURE 9A

corr_20days = srisk_analysis_data %>%
  filter(series %in% c("aggregate_srisk_truncated_bank", "aggregate_srisk_truncated_nonbank")) %>%
  create_correlations(bank_index = "aggregate_srisk_truncated_bank", nonbank_index = "aggregate_srisk_truncated_nonbank", rollcorr_param = 20, lag_param = 0) %>%
  date_to_period_mapping() %>%
  group_by(nonbank, Period, Window, period_number) %>% 
  summarize(median = median(cor, na.rm = T) %>% round(3)) %>%
  ungroup() %>% 
  arrange(period_number) %>%
  select(Period, Window, median) %>%
  setNames(c("period", "window", "median_correlation")) 
corr_20days

fig_9a <- 
  ggplot(corr_20days, aes( x=period, y = median_correlation)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(x = "", y = "", title = "20-day Rolling Correlations between Bank and NBFI % Changes in Sector-wide SRISK") +
  coord_cartesian(ylim = c(0.55, 0.95)) 

fig_9a


# FIGURE 9B

# Join with s&p 500 returns and recover abnormal returns
granger_analysis_data = srisk_analysis_data %>%
  filter(series %in% c("equity_returns_bank_equal_weighted", "equity_returns_nonbank_equal_weighted")) %>%
  inner_join(sp500, by = "date") %>% extract_alpha() %>% 
  filter(! is.na(abnormal_return)) %>%
  select(date, series, abnormal_return) %>%
  pivot_wider(names_from = "series", values_from = "abnormal_return")

# Loop through every rolling Granger causality test
granger_analysis_data$bank_cause_nonbank = NA; granger_analysis_data$nonbank_cause_bank = NA
for (i in 91:NROW(data)) {
  subset = granger_analysis_data %>% slice((i-90):i)
  granger_analysis_data$bank_cause_nonbank[i] = granger_function(subset$equity_returns_bank_equal_weighted, subset$equity_returns_nonbank_equal_weighted, max_lags = 10, optimal_lags = T)
  granger_analysis_data$nonbank_cause_bank[i] = granger_function(subset$equity_returns_nonbank_equal_weighted, subset$equity_returns_bank_equal_weighted, max_lags = 10, optimal_lags = T)
  print(i)
}

fig_9b <- granger_analysis_data %>% mutate(significant_bank_cause_nonbank = ifelse(bank_cause_nonbank < 0.1, 1, 0),
                significant_nonbank_cause_bank = ifelse(nonbank_cause_bank < 0.1, 1, 0)) %>%
  date_to_period_mapping() %>%
  group_by(Period) %>%
  summarize(bank_cause_nonbank = mean(significant_bank_cause_nonbank, na.rm = T),
            nonbank_cause_bank = mean(significant_nonbank_cause_bank, na.rm = T)) %>%
  ungroup() %>%
  mutate(bank_cause_nonbank = paste0((round(bank_cause_nonbank, 2))*100, "%"),
         nonbank_cause_bank = paste0((round(nonbank_cause_bank, 2))*100, "%")) %>%
  setNames(c("Period", "Bank Causes Nonbank", "Nonbank Causes Bank")) %>%
  gt() %>%
  tab_spanner(label = "90 Day Window", columns = c("Bank Causes Nonbank", "Nonbank Causes Bank"))

fig_9b